using System;

namespace BOL.Maths.LinearAlgebra
{
	/// <summary>
    /// Cholesky Decomposition.
	/// For a symmetric, positive definite matrix A, the Cholesky decomposition
	/// is an lower triangular matrix L so that A = L*L'.
	/// If the matrix is not symmetric or positive definite, the constructor
	/// returns a partial decomposition and sets an internal flag that may
	/// be queried by the IsSPD() method.
	/// </summary>
	public sealed class CholeskyDecomposition : DecompositionBase
	{
		#region Private variables

        private double[][] _l;
        private bool _isSymmetric;
        private bool _isPositiveDefinite;

		#endregion

		#region Public Properties
		
        /// <summary>Returns triangular factor.</summary>
        public double[][] L { get { return _l; } }

        /// <summary>Returns whether the matrix is symmetric.</summary>
        public bool IsSymmetric { get { return _isSymmetric; } }

        /// <summary>Returns whether the matrix is positive definite.</summary>
        public bool IsPositiveDefinite { get { return _isPositiveDefinite; } }

        #endregion

        #region Constructor

        /// <summary>Cholesky algorithm for symmetric and positive definite matrix.</summary>
		public CholeskyDecomposition(double[][] matrix) : base(matrix)
		{
            if (!MatrixDouble.IsSquare(matrix))
                throw new ArgumentException("Matrix must be square.");

            _l = new double[_m][];
            for (int i = 0; i < _m; i++)
                _l[i] = new double[_m];
            _isSymmetric = true;
            _isPositiveDefinite = true;

            InitilizeL(matrix);
		}

        private void InitilizeL(double[][] a)
        {
            double[] row_j, row_k;
            double d, s;

            for (int j = 0; j < _m; j++)
            {
                row_j = _l[j];
                d = 0.0;

                for (int k = 0; k < j; k++)
                {
                    row_k = _l[k];
                    s = 0.0;

                    for (int i = 0; i < k; i++)
                        s += row_k[i] * row_j[i];

                    row_j[k] = s = (a[j][k] - s) / _l[k][k];
                    d = d + s * s;
                    _isSymmetric &= (a[k][j] == a[j][k]);
                }

                d = a[j][j] - d;
                
                _isPositiveDefinite &= (d > 0.0);
                _l[j][j] = Math.Sqrt(Math.Max(d, 0.0));
                for (int k = j + 1; k < _m; k++)
                    _l[j][k] = 0.0;
            }
        }
		
		#endregion

		#region Public Methods
		
		/// <summary>
        /// Solve A * X = B
		/// <param name="b">A Matrix with as many rows as A and any number of columns.</param>
		/// <returns>X so that L * L' * X = B</returns>
		public double[][] Solve(double[][] b)
		{
            if (b.Length != _m)
                throw new ArgumentException("Matrix row dimensions must agree.");

			if (!_isSymmetric)
				throw new SystemException("Matrix is not symmetric.");

			if (!_isPositiveDefinite)
				throw new SystemException("Matrix is not positive definite.");
			
            int i, j, k;
            int mb = b.Length;
            int nb = b[0].Length;
			// copy right hand side.
			double[][] x = new double[mb][];
            for (i = 0; i < mb; i++)
                Array.Copy(b[i], x[i], nb);
			
			// solve L * Y = B;
			for (k = 0; k < _m; k++)
			{
				for (i = k + 1; i < _m; i++)
					for (j = 0; j < mb; j++)
						x[i][j] -= x[k][j] * _l[i][k];

				for (j = 0; j < mb; j++)
					x[k][j] /= _l[k][k];
			}
			
			// solve L' * X = Y;
			for (k = _m - 1; k >= 0; k--)
			{
				for (j = 0; j < mb; j++)
					x[k][j] /= _l[k][k];

				for (i = 0; i < k; i++)
					for (j = 0; j < mb; j++)
						x[i][j] -= x[k][j] * _l[k][i];
			}

            return x;
		}

		#endregion
	}
}
